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Abstract When granular systems are modeled by fric- 
tionless hard spheres, particle-particle collisions are con- 
sidered as instantaneous events. This implies that while 
the velocities change according to the collision rule, the 
positions of the particles are the same before and after 
such an event. We show that depending on the material 
and system parameters, this assumption may fail. For 
the case of viscoelastic particles we present a univer- 
sal condition which allows to assess whether the hard- 
sphere modeling and, thus, event-driven Molecular Dy- 
namics simulations are justified. 

Keywords Granular Gases • Hard Sphere Model • 
Coefficient of Normal Restitution • Viscoelastic 
Spheres • event-driven Molecular Dynamics 
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1 Introduction 

Hard sphere modeling of granular systems assumes that 
the dynamics of the system may be described as a se- 
quence of instantaneous events of binary collisions. In 
between the collisions the particles move freely along 
straight lines, or ballistic trajectories in presence of ex- 
ternal fields like gravity. The hard-sphere model of par- 
ticle collisions is the foundation of both Kinetic The- 
ory of granular matter, based on the Boltzmann equa- 
tion, e.g. [5 ,3 ,14 , and event-driven Molecular Dynam- 
ics (eMD) of granular matter, e.g. 
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In hard sphere approximation, the inelastic collision 
of frictionless spheres i and j located at r^ and r j travel- 
ing at velocities r^ and Tj is, thus, characterized by the 
collision rule describing the instantaneous exchange of 
momentum between the colliders, 
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with the unit vector e ? 
index describes values just before the collision, primed 
values describe postcollisional values. The inelasticity is 
characterized by the coefficient of (normal) restitution 
e. 

The instantaneous character of the collisions implies 
that as the result of a collision only the velocities of the 
particles change but not their positions, r- = r'- = rf- 
and, thus, e' r = With this, Eq. ([!]) turns into 

ft - r;.) • e° = -s HS (r? - r°) • e° (2) 

which allows to compute the postcollisional velocities 
successively for all collisions in the system, which is the 
basic idea of event-driven Molecular Dynamics (EMD). 
Provided the system may be described as hard spheres 
undergoing instantaneous collisions, EMD may be by 
orders of magnitude more efficient than ordinary MD 
integrating Newton's equation of motion. Recently, ex- 
tremely efficient algorithms for EMD simulations have 
been developed, e.g. pQ. 

As physical particles are not perfectly hard but the 
collision is governed by finite interaction forces, the 
hard sphere model is an idealization whose justifica- 
tion may be challenged. Especially in view of its im- 
portance for Kinetic Theory and numerical simulation 
techniques. In particular, for finite duration of the col- 
lisions the unit vector e r may rotate during a collision 
by the angle 



a 



arccos (e 
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invalidating Eq. ([2| and, therefore, the hard-sphere ap- 
proximation. While this angle is negligible for approx- 
imately central impacts of relatively stiff spheres, it is 
not for oblique impacts of soft spheres [2]. 

Within this work we quantify under which condi- 
tions and to what extend the condition, e' r « of the 
hard sphere assumption fails. Aim of the present pa- 
per is to provide a universal condition which allows for 
arbitrary collisions of arbitrary elastic spheres to assess 
whether the hard sphere model is acceptable for the 
description of particle collisions. Thus, we discriminate 
wether the hard sphere model is acceptable for systems, 
characterized by (i) a set of material parameters, (ii) 
particle sizes and (iii) a typical (thermal) impact ve- 
locity. To generalize our result to the case of inelastic 
collisions, we show that regarding the rotation angle a, 
elastic spheres are the marginal case, that is, if e' r ~ 
holds true for elastic particles, it certainly holds true 
for inelastic particles. 



2 Collision of spheres 

Consider two colliding spheres with the masses mi and 
rrij located at and rj and traveling with velocities 
and r.. With the interaction force F, their motion is 



described by 

m e frr = F , MR = 



where 



R 



rriiXi 



^eff 



rn.rnj 



(4) 



(5) 



are the center of mass coordinate, the relative coordi- 
nate and the effective mass, respectively. The center of 
mass moves due to external forces such as gravity and 
separates from the relative motion which in turn con- 
tains the entire collision dynamics. 

For frictionless particles, the interaction force ex- 
clusively acts in the direction of the inter-center unit 
vector, F = F n e ri that is, there is no tangential force 
and, thus, the particles' rotation is not affected by the 
collision. During the collision the (orbital) angular mo- 
mentum is conserved which allows for the definition of 
a constant unit vector: 



L = m e ff r x r = Le-L . 

Thus, with the coordinate system spanned by 
s — 



(6) 



(7) 



and with its origin in the center of mass R, the collision 
takes place in the e^-e^-planeQln the collision plane we 

1 For central collisions we have L = 0. In this case e z may 
be any unit vector perpendicular to e x , (e x • e z = 0). 




Fig. 1 Illustration of the used polar coordinates (see text) 



formulate the equation of motion in polar coordinates 
{r, cp} (see Fig. [I]): 

m eff rV = ^, mefff = F c + F n = m e ^r(p 2 + F n , (8) 

with the centrifugal force F c . Together with the inital 
conditions 



r(0) = r°, r(0) = r°, <p(0) = , 



(9) 



Eq. (|8| fully describes the collision dynamics for an 
arbitrary normal force F n . The collision terminates at 
time t = r where [T8lfT9] 



r(r) > and F n = 0. 



(10) 



Inserting the first equation of Eq. (|8| into the second, 
we obtain 



i? 

m e ffr = ^ + F n 



(11) 



which fully governs the radial dynamics of the problem. 

Note that in contrast to earlier work [18,19 where 
the coefficient of normal restitution was derived from 
force laws F n here we allow the normal vector e r to ro- 
tate during the collision and do not neglect the resulting 
centrifugal force. 

Since for any finite interaction forces, the duration 
r of a collision is finite, for non-central collisions, L ^ 
0, during the collision the spheres rotate around their 
center of mass, that is, e' r ^ e® and see Eq. (|3|. 

It is frequently stated that the hard sphere approx- 
imation and thus event-driven simulations are always 
justified for dilute systems where the mean free flight 
time of the particles is large compared to the typical col- 
lision time. Obviously, this condition is insufficient. It 
may be shown that the characteristics of dilute granular 
gases such as the coefficient of self diffusion sensitively 
depends on the rotation of the unit vector [T2] . 
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3 Elastic Spheres 

3.1 Dimensonless Equation of Motion 

The collision of elastic spheres obeys Hertz' contact 
force [6], 



p{l-rf'\ I 



Ri + Rj 



(12) 



where I denotes the distance between the particle cen- 
ters at the moment of impact. The quantity £ = I — r 
is often referred to as the deformation or mutual com- 
pression. The elastic constant p reads 



(13) 



3(1 -is 2 ) 



where Y, v and R e ^ stand for the Young modulus, the 
Poisson ratio and the effective radius R e ^ = RiRj /(Ri + 
Rj), respectively. 

Writing the general equation of motion Eq. ^ with 



the force given by Eq. (12) and measuring length in 



units of X and time in units of T [19] 
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According to Eqs. ( 15]) and ([16]) together with the initial 
conditions Eq. dl9), the binary collision of frictionless 
elastic spheres is described by only two free parameters: 
I and c cn . 



3.2 Rotation of the Normal Vector 



T L 



^ ~ X 2 m eff * 

The scaled initial conditions read 
cp(0) = 0, f (0) = I and r(0) = 1. 



We solve the equations of motion (15), (16) and (19) 



to obtain the rotation a of the unit vector, e r , given 
by Eq. ([3| as a consequence of the collision of elastic 
spheres. This rotation occurs for oblique impacts and is 
commonly neglected in hard sphere simulations as well 
as in the Kinetic Theory of granular gases. Obviously, 



a depends on the material properties, the particle sizes 
and on the geometry of the collision as sketched in Fig. 
[2j The condition a w will be interpreted as an index 
for the justification of the hard-sphere approximation 
for a given system. 




Fig. 2 Eccentric binary collision of spheres. The sketched 
situation corresponds to the eccentricity d/l ~ 0.8 where the 
rotation angle drawn in Fig. [3] adopts its maximum. 



To illustrate the fact that the rotation of the unit 
vector is by far not a small effect even for rather com- 
mon systems, in Fig. [3] we plot the angle a as a func- 
tion of the impact eccentricity d/l (see Fig. [2J. The 
system parameters (in physical units) are: radii R\ = 
i?2 =0.1 m, material density p m = 1140kg/m 3 , Young 
modulus Y = 10 T N/m 2 , Poisson ratio v = 0.4, impact 
velocity 20m/s. 




Fig. 3 Rotation angle a of the unit vector e r as a function 
of the impact eccentricity d/l (see Fig. [2]) for rubber spheres, 
parameters specified in the text. The marked area shows the 
interval where a > 30° corresponding to 65% of all collisions 
when molecular chaos is assumed. 



As expected, the rotation vanishes for central col- 
lisions. The rotation adopts its maximum for d/l « 
0.8 (this situation corresponds to the sketch in Fig. |2| 
where it can easily reach values of a w 40°. The po- 
sition of the maximum may surprise since in the Ki- 
netic Theory it is frequently assumed that if at all only 
rare glancing collisions might deserve a special consid- 



4 



Patric Miiller, Thorsten Poschel 



eration. Assuming molecular chaos, that is, e = d/l is 
distributed as dp(e) = 2e de, and the parameters given 
above, about 65% of the collisions lead to a rotation an- 
gle a > 30° (marked interval in Fig. [3|. Consequently, 
the rotation of the unit vector e r is a very significant 
effect for granular gases. 



3.3 Universal Description of the Rotation Angle 

For elastic particles the dimensionless equation of mo- 
tion of the collision (Eqs. (p]), ([16]) and (l9|» is fully 



specified by two independent parameters, I and c^, de- 



fined in Eqs. (17) and (18). Therefore all material and 



system parameters may be mapped to a point in the 
(I , c^,)-space. 

The rotation angle a can be determined by the fol- 
lowing procedure: 

1. Determine the dimensionless parameters: 
{Y, v, R, p m , v, d/l} {I, c^} 

2. Solve numerically the equations of motion, 

Eqs. ( 15|16|19 ) for < t < r where r is the time 
when the collision terminates, r is determined by 
the conditions f(r) = and f(r) > (see Eq. ([Toj)) . 

3. The rotation angle is obtained from a = ip(r). 

We performed this procedure for a wide range of rel- 
evant (physical) parameters given in Table [l] In dimen- 





unit 


min. 


max. 




Y 


[10 9 N/m 2 ] 


0.01 


100 


Young's modulus 


V 




0.2 


0.5 


Poisson ratio 


R 


[m] 


0.001 


0.1 


particle radius 


Pm 


[kg/m 3 ] 


250 


3250 


material density 


V 


[m/s] 


0.001 


25 


impact velocity 


d/l 




0.01 


0.99 


eccentricity 



Table 1 Physical parameter space scanned to obtain Fig. [4j 
For the definition of the impact velocity and the eccentricity 
we refer to Fig. [2] 



sionless variables, this range corresponds to the interval 



2.12 < I < 1.8 • 10 9 , 2.12 • 10" 2 < c<p < 1.26 • 10 9 . (20) 

Figure ^ shows the rotation angle a as a function of I 

and on a double logarithmic scale. 

Using the definitions of f, and L, Eqs. ( 17|l8|6 ) 
and 
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Fig. 4 Rotation angle a as function of I and . Grey regions 
indicate points which do not correspond to any combination 
of parameters given in Tab. [l] Only combinations inside the 
dashed region may lead to a noticeable rotation angle. 



which follows from the definitions, Eq. (14), and geom- 
etry (see Fig. one obtains 

2 

- 1 



In/ 



1 



ln(c^) + ^ln 



(22) 



This equation provides some insight into the structure 
of Fig. [4] and allows for a more intuitive presentation 
of the result. For fixed eccentricity d/l due to Eq. p2| ), 
in the double logarithmic scale used in Fig. [ij I is a 
linear function of with slope 1. That is, all collisions 
taking place at the same impact eccentricity d/l are 
located on a straight line of slope 1 in the (lnc^,ln/)- 
space. The position along this line is then determined 
by the remaining system parameters. 

The chosen interval, 0.01 < d/l < 0.99, see Tab. [l] 
implies that the intercept, of all possible straight lines 



given by Eq. (22) is bound to the range 
-1.95 < ln[ < 4.61, 



(23) 



which explains the stripe structure of the data in Fig. |1J 
All (In Z, In c^)-pairs outside the colored stripe cannot be 
adopted for any combination of the parameters listed 
in Tab. [l] which is indicated by the gray areas in Fig. [4] 
Figure [4] indicates that among all studied combina- 
tions of parameters only those for — 3 ^ lnc^, ^ 9 and 
1 ;$ In I ;$ 8 (dashed region in Fig. [4jTmay lead to a 
noticeable rotation angle a or a significant deviation 
from the hard sphere model respectively. Therefore, we 
study this region with a higher resolution of and /, 
see Fig. [5] In order to avoid drawing the irrelevant gray 
regions, we plot the data over In — In c™ instead of 



lnc^ with 



lnc™ = lnZ 



In 



(f) 



2 

min 



In [-4.61. (24) 
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as obtained from Eq. (22) with (d/Z) m in taken from Tab. 
[I] (see illustration of lnc^ — lnc™ in Fig. [I]). 



We specify a collision by 




40 g> 



Fig. 5 Rotation angle a as function of I and c^/c™ 111 . The 
figure shows the region of the parameter space where a adopts 
noticeable values (dashed region in Fig.^J. Additionally sev- 
eral isolines of constant rotation angle a are shown. 



The isolines of constant rotation angle a drawn in 
Fig. [5] indicate that there is a rather sharp transition 
between the regions where a « and a ^> in the 
(lnc^, ln/)-space. Hence, regarding the rotation of the 
unit vector e r , the regions in the parameter space where 
the hard sphere model is a justifiable approximation 
are clearly separated from those, where the hard sphere 
approximation is questionable. 



3.4 Confidence Regions of the Hard Sphere Model 

For practical applications one might wish to know whether 
a given set of material and system parameters allows 
for a hard-sphere description. Besides other criteria, 
the maximum rotation angle which can be expected for 
these parameters, is an important criterion. For the fol- 
lowing we assume that the rotation angle a c is marginally 
acceptable for the hard sphere approximation and pro- 
vide a simple approximate method do decide whether 
the given system fulfills this criterion. 

Fig. [5] shows that for rotation angles up to about 
15°, the isolines of constant rotation angle are approxi- 
mately straight lines of slope m ~ 0.84 on average. The 
corresponding intercept t OLc increases with the isoline 
value a c . From Fig. [5] we obtain t\o « 0.9, £50 w —0.29, 
tio° ~ -0.74 and t lb o « -0.85. 



In 



= 4.61 



In 



3(1 - v 2 )m e ff_ 



2/5 



vy/i - (d/iy 



-4/5 



and define 

D ac =\nl-(m\n(c^/c™ n )+t ac ) . 



(25) 



(26) 



D 1 > indicates that the maximally expected rotation 
angle is smaller than a c , that is, the hard sphere model 
is acceptable for this situation. 



4 Inelastic Spheres 

4.1 Equation of Motion 

The main conclusion of this Section will be that in- 
elastic interaction forces which are, perhaps, the most 
characteristic feature of granular materials, do not lead 
to an increase of the rotation angle a as compared with 
the elastic case detailed in the previous Section. Here, 
we exemplary discuss a particular dissipation mecha- 
nism, the viscoelastic model which is widely used for 
modeling granular systems, e.g. [7 , 20lfT7] . Many other 
dissipative interaction forces as plastic deformation, lin- 
ear dashpot damping, etc. lead to very similar results. 

The collision of viscoelastic spheres is characterized 
by the interaction force [4] 

p(l - r) 3/2 - -AprVT^V (27) 



ridis 



with the dissipative constant A being a function of 
the elastic and viscous material parameters [4 and the 
other parameters as described before. The collision ter- 
minates at time r when r(r) < and f(r) = 0, corre- 
sponding to purely repulsive interaction, see [19]. The 
dissipative part, i^ ls , was first motivated in [8] and 
then rigorously derived in [4] and [IT , where only the 
approach in [4 leads to an analytic expression of the 
material parameter A. 



We apply the same scaling as in Sec. 3.1 to obtain 



d(f c 
dt ~ r 
_ „ 
dt 2 ~ T \dt 



dip 



+ I 



3/2 



Cdit 



r — 
dt 



(28) 



with the definitions and initial conditions given in Eqs. 
([1T| [18) [19]) and additionally 

3 pA 
2 m e ff ' 



c dis jVXT ; 
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(29) 
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In contrast to the case of elastic spheres discussed in 
Section [3] for inelastic frictionless spheres we need three 
independent parameters to describe their collisions, Z, 
c^ and c dis . 



4.2 The Role of Inelasticity 

To study the dependence of the rotation angle a on 
the inelasticity, we repeat the computation shown in 
Sec. |3.2| (same elastic parameters) for inelastic collisions 
where A ^ 0. Figure [6] shows the rotation of the unit 
vector e r during inelastic collisions over the eccentricity 
d/l (see Fig. |2| for various dissipative constants A. 




Fig. 6 Rotation angle a of the unit vector e r as function 
of the eccentricity for various dissipative constants A. The 
elastic parameters are the same as for Fig. [3] 



In Fig. [7] we additionally fix d/l = 0.5 to plot the ro- 
tation a as a function of the dissipative constant A. To 
provide a more vivid quantity for the inelasticity of the 
collision, we give a also as a function of the coefficient of 
restitution e corresponding to a central collision at the 
chosen impact velocity v n = y/3/2 • 20 m/s « 17.3 m/s 

Disregarding for a moment the centrifugal force, 
the dependence of the moment of inertia on l(t) and 
the fact that the final deformation Z(r) depends on 
A [15] , the decreasing function a (A) or a(e) may be 
understood essentially from the fact that the duration 
of the contact is a decreasing function of inelasticity, 
dr(A)/dA < 0. Thus, the smaller the coefficient of resti- 
tution the shorter lasts the contact and the smaller is 
the rotation angle during contact. This explanation is 
certainly oversimplified and serves only as a motivation 
to understand qualitatively the behavior of a(A). 

As shown qualitatively in Fig. [6] and quantitatively 
in Fig. [7] for all eccentricities the rotation angle adopts 




0.002 0.004 0.006 

A (1/s) 



0.008 



0.01 



Fig. 7 Rotation angle a over the dissipative parameter A 
(lower scale, full line) and over the coefficient of normal resti- 
tution e (upper scale, dashed line) for d/l = 0.5. 



its maximum for A = 0, corresponding to elastic colli- 
sions, £ = 1. 



5 Conclusion 

For all real materials the collision of particles implies 
a rotation of the inter-particle unit vector e r during 
the time of contact r by a certain angle a. This rota- 
tion is neglected in Kinetic Theory of granular systems 
as well as in event-driven Molecular Dynamics simula- 
tions relying both on the hard-sphere model of granu- 
lar particles. Therefore, to justify the application of the 
hard-sphere model, one has to assure that the rotation 
angle is negligible for the given system parameters. In 
the present paper, we reduce the problem of oblique 
elastic collisions to two independent parameters, I and 
and compute the rotation angle a as a function 
of these parameters. The result is universal, that is, a 
is known for any combination of material parameters 
(Young modulus Y, Poisson ratio /i, material density 
p m ) and system parameters (particle radii R, impact 
eccentricity d/l, and impact velocity v). 

For dissipative collisions characterized by the coef- 
ficient of restitution, < e < 1, we show that the ro- 
tation angle is smaller than for the corresponding elas- 
tic case where all parameters are the same, except for 
s = 1. Therefore, to assess whether the rotation angle is 
small enough to justify the hard sphere approximation 
for a given system of dissipative particles, it is sufficient 
to consider the corresponding system of elastic particles 
discussed in Sec. 03 

For convenient use of our result we provide a univer- 
sal lookup table and the corresponding access functions 
(see Online Resource [TO]). The angle of rotation for a 
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given situation can be either obtained using the dimen- 18. 



sionless variables, a(c^, I), obtained from Eqs. (17) and 



(18) together with Eqs. (13) and (14) for the present 



material and system parameters, or by providing the 
physical parameters directly. 

Concluding we consider our result as a tool to assess 
whether the Kinetic Theory description of a granular 
system on the basis of the Boltzmann equation and/or 
its simulation by means of highly efficient event-driven 
Molecular Dynamics is justified. 



19. 



20. 



Schwager, T., Poschel, T.: Coefficient of restitution and 
linear-dashpot model revisited. Granular Matter 9, 465 
(2007) 

Schwager, T., Poschel, T.: Coefficient of restitution for 
viscoelastic spheres: The effect of delayed recovery. Phys. 
Rev. E 78 (2008) 

Stevens, A.B., Hrenya, CM.: Comparison of soft-sphere 
models to measurements of collision properties during 
normal impacts. Powder Technology 154, 99 (2005) 



Acknowledgements The authors gratefully acknowledge the 
support of the Cluster of Excellence 'Engineering of Advanced 
Materials' at the University of Erlangen-Nuremberg, which is 
funded by the German Research Foundation (DFG) within 
the framework of its 'Excellence Initiative'. 



References 

1. Bannerman, M.N., Sargant, R., Lue, L.: An o(n) general 
event-driven simulator: DYNAMO. J. Comp. Chem. (in 
press) (2011) 

2. Becker, V., Schwager, T., Poschel, T.: Coefficient of tan- 
gential restitution for the linear dashpot model. Phys. 
Rev. E 77, 011304 (2008) 

3. Brilliantov, N.V., Poschel, T.: Kinetic Theory of Gran- 
ular Gases. Oxford graduate texts. Oxford University 
Press, Oxford (2010) 

4. Brilliantov, N.V., Spahn, F., Hertzsch, J.M., Poschel, T.: 
Model for collisions in granular gases. Phys. Rev. E 53, 
5382 (1996) 

5. Goldhirsch, I.: Rapid granular flow. Ann. Rev. Fluid 
Mech. 35, 267 (2003) 

6. Hertz, H.: iiber die Beriihrung fester elastischer Korper. 
Journal fur reine und angewandte Mathematik 92, 156 
(1881) 

7. Kruggel-Emden H. ans Simek, E., Rickelt, S., Wirtz, S., 
Scherer, V.: Review and extension of normal force mod- 
elsfor the discrete element method. Powder Technology 
171, 157 (2007) 

8. Kuwabara, G., Kono, K.: Restitution coefficient in a col- 
lision between two spheres. Jpn. J. Appl. Phys 26, 1230 
(1987) 

9. Lubachevsky, B.D.: How to simulate billards and similar 
systems. J. Comp. Phys. 94(2), 255 (1991) 

10. TODO: Please insert link to Online Resource. 

11. Morgado, W.A.M., Oppenheim, I.: Energy dissipation for 
quasielastic granular particle collisions. Phys. Rev. E 55, 
1940 (1997) 

12. Miiller, P., Poschel, T.: in preparation 

13. Miiller, P., Poschel, T.: Collision of viscoelastic spheres: 
Compact expressions for the coefficient of restitution. 
Phys. Rev. E (in press) (2011) 

14. Poschel, T., Brilliantov, N. (eds.): Granular Gas Dynam- 
ics. Lecture Notes in Physics. Springer, Berlin (2003) 

15. Poschel, T., Schwager, T.: Computational Granular Dy- 
namics. Springer, Berlin (2005) 

16. Rapaport, D.C.: The event scheduling problem in molec- 
ular dynamic simulation. Journal of Computational 
Physics 34, 184 (1980) 

17. Schafer, J., Dippel, S., Wolf, D. E.: Force schemes in sim- 
ulations of granular materials. J. Phys. I France 6, 5 
(1996) 



